/*
This is a very brief program in order to test the ewald code.

In the header file 'ewald.h' there is a maximal limit for the
reciprocal space cutoff, which (if required) has to be adjusted.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fstream>
#include <iostream>
#include <iomanip>

using namespace std;

#define maxkNew 100


/* maximal precision (used in cutoff) */
#define PREC 1.0e-30

/* precision of alpha */
#define ALPHA_OPT_PREC 1.0e-10

/* using the formula of Benson-Mackenzie on gets for the Madelung
constant of the NaCl lattice: */
const double madelung_NaCl = -1.747564594633182190636212035544397403481;

/* Berechnet das Quadrat von double x, gibt double zurueck. */
double SQR(double x) { return x*x; }

/* Wert von double x runden, gibt double zurueck. */
double dround(double x) { return floor(x+0.5); }

#define     PI     3.14159265358979323846264338328
#define   wupi     1.77245385090551602729816748334

/* constant to use when alpha is to be optimized */
#define ALPHA_OPT 0.0

#include "utils.h"
#include "directo.h"
#include "ewdClass.h"

int main(void)
{
	int i;
	/*Numero de particulas*/
	//int npart = 512;
	int npart =10;
	/*Posiciones de las particulas*/
	double **pos;
	/*Arrays de fuerzas en espacio directo y reciproco*/
	double **forceR, **forceK, **forceUsual, **forceUsualReplica;
	/*Lado del cubo de simulacion*/
	double lbox = 10.0;
	double s1 = 0.50*lbox;
	/*Cutoff del potencial en espacio directo (tiene que ser menor que 0.5*lbox)*/
    double rcut  = 4.9;
	//double rcut  = 0.499*lbox;
	double rcut2;
	/*Cutoff del potencial en el espacion reciproco*/
	int    kmax  = 10;
	/*Coeficiente de corte de Ewald*/
	double Ewaldalpha=7.7;
	/*Suma de la carga cuadratica (Coulomb)*/
	double Q2 = 8.0;
	/*Array de carga de partiMiArchivo>>pos[i][0]>>pos[i][1]>>pos[i][2];culas (Coulomb)*/
	double *Q;
	/*Prefactor del potencial Coulombiano (Coulomb) */
	double coulomb_prefactor = 1.0;
	/*Constante dielectrica del entorno*/
	double epsilon = 1.0;   /* = Vacio */
	/* real space cutoff. Note that this must be smaller than L/2, if
	minimum image convention is to be used! */
	/*Energia del sistema*/
	double energia = 0.0;
	double energiaDirecta=0.0;
	double energiaDirectaReplica=0.0;
	void loadconfig(double **pos);
	void loadconfig(double **pos, int npart,double lbox);
	double eDirecta, eReciproca;

for (npart=5;npart<20;npart++){
	/*Creamos los arrays*/
	pos = new double*[npart];
	forceR=new double*[npart];
	forceK=new double*[npart];
	forceUsual=new double*[npart];
	forceUsualReplica=new double*[npart];
	Q = new double[npart];
	for(int i=0;i<npart;i++){
		pos[i]=new double[3];
		forceR[i]=new double[3];
		forceK[i]=new double[3];
		forceUsual[i]=new double[3];
		forceUsualReplica[i]=new double[3];
	}



	/*Creamos el objeto de suma de Ewald*/
	ewald *ewd=new ewald();

	/*Posiciones de las particulas y cargas. Corresponden a una red NaCl*/
//	pos[0][0] =  0;  pos[0][1] =  0;  pos[0][2] =  0;  Q[0] =  1.0;
//	pos[1][0] = s1;  pos[1][1] = s1;  pos[1][2] =  0;  Q[1] =  1.0;
//	pos[2][0] = s1;  pos[2][1] =  0;  pos[2][2] = s1;  Q[2] =  1.0;
//	pos[3][0] =  0;  pos[3][1] = s1;  pos[3][2] = s1;  Q[3] =  1.0;
//
//	pos[4][0] = s1;  pos[4][1] =  0;  pos[4][2] =  0;  Q[4] =  -1.0;
//	pos[5][0] =  0;  pos[5][1] = s1;  pos[5][2] =  0;  Q[5] =  -1.0;
//	pos[6][0] =  0;  pos[6][1] =  0;  pos[6][2] = s1;  Q[6] =  -1.0;
//	pos[7][0] = s1;  pos[7][1] = s1;  pos[7][2] = s1;  Q[7] =  -1.0;


	rcut2 = SQR(rcut);
    loadconfig(pos,npart,lbox);
	/* initialize the Ewald-algorithm */
	//Ewaldalpha=ALPHA_OPT;

	/*Importamos el archivo de configuracion*/
	//loadconfig(pos,&lbox);

	/*Inicializacion para el calculo de potencial dipolar*/
	ewd->EwaldDipolar_init(lbox,npart,pos, forceR, forceK, Ewaldalpha,rcut,kmax,epsilon);
//	ewd->EwaldCoulomb_init(lbox, npart, pos, forceR, forceK,Q,  Ewaldalpha, rcut,
//		kmax,  coulomb_prefactor, Q2, epsilon);

	/*Inicializamos las fuerzas*/
	for (i = 0; i < npart; i++){
		for(int k=0;k<3;k++){
			forceR[i][k]= 0.0;
			forceK[i][k]= 0.0;
		}
	}

	/* Calculamos la contribucion en el espacio real de las energias/fuerzas dipolares*/
	eDirecta = ewd->EwaldDipolar_r_space(1);
	/* Calculamos la contribucion en el espacio real de las energias/fuerzas Coulombianas*/
	//eDirecta = ewd->EwaldCoulomb_r_space(1);

	/* Calculamos la contribucion en el espacio reciproco de las energias/fuerzas dipolares*/
	eReciproca = ewd->EwaldDipolar_k_space(1);
	/* Calculamos la contribucion en el espacio reciproco de las energias/fuerzas Coulombianas*/
	//eReciproca = ewd->EwaldCoulomb_k_space(1);
	/*Calcula la energia y fuerza dipolar de forma usual*/
	energiaDirecta=dipolarEnergyForces(pos,ewd->lbox,forceUsual,ewd->npart);
	//energiaDirecta=CoulombEnergyForcesReplication(pos,Q,ewd->lbox,ewd->npart);

	/*Calculamos la energia y fuerza con replicas*/
	energiaDirectaReplica= dipolarEnergyForcesReplication(pos, lbox, forceUsualReplica, npart);

//    printf("Energia Directa: %6.5f\n",eDirecta);
//    printf("Energia Reciproca: %6.5f\n",eReciproca);
//	printf("Energia Ewald: %6.5f\n",eDirecta+eReciproca);
//	printf("Energia con replicas: %6.5f\n",energiaDirectaReplica);
//	printf("Energia Directa: %6.5f\n",energiaDirecta);
	printf("%d  %f  %f  %f  %f  \n",npart,eDirecta,eReciproca,(eDirecta+eReciproca),energiaDirectaReplica);

	/*Eliminamos los objetos creados*/
	delete ewd;
	for(int i=0;i<3;i++){
		delete pos[i];
		delete forceR[i];
		delete forceK[i];
		delete forceUsual[i];
		delete forceUsualReplica[i];
	}
	delete []pos;
	delete []forceR;
	delete []forceK;
	delete []forceUsual;
	delete []forceUsualReplica;
	delete Q;
}
	return 0;
}

void loadconfig(double **pos)
{
    ifstream MiArchivo ("/home/jonk/pos.dat",ios::in);
	char line[150];
	char dummy[20];
	int npart=100;
	if (MiArchivo.is_open())
	{
	    for(int i=0;i<npart;i++)
        {
            MiArchivo>>pos[i][0]>>pos[i][1]>>pos[i][2];

        }
        MiArchivo.close();
	}

}

void loadconfig(double **pos, int npart,double lbox)
{
	for(int i=0;i<npart;i++)
	{
		for(int k=0;k<3;k++)
		{
			//pos[i][k]=lbox*(double(rand())/RAND_MAX-0.5);
			pos[i][k]=(double(rand())/RAND_MAX-0.5);
		}
	}
}







